package fem2.material;

import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.math.LinAlg;

import java.util.ArrayList;
import java.util.Map;
import java.util.TreeMap;

import math2.FunctionRnToR;
import math2.LinAlgUtilities;
import math2.MutableDouble;
import fem2.Debugger;
import fem2.Element;
import fem2.element.StructuralElement;
import fem2.enu.EchoLevelType;
import fem2.enu.State;
import fem2.enu.VariableName;
import fem2.mesh_and_topology.IntegrationPointIndexer;
import fem2.mesh_and_topology.SpatialCriteria;
import fem2.mesh_and_topology.SphereSpatialCriteria;

/**
 * implicit/explicit version of isotropic damage material model
 * 
 * @author hbui
 * 
 */
public class IsotropicDamageMMimplexNonlocal extends IsotropicDamageMMimplex {

	private FunctionRnToR weightFunction;
	private double interactionRadius;

	private IntegrationPointIndexer indexer;
	private Map<Element, Map<Integer, Map<Element, ArrayList<Integer>>>> integrationPointsMapping = new TreeMap<Element, Map<Integer, Map<Element, ArrayList<Integer>>>>();

	public IsotropicDamageMMimplexNonlocal(double E, double nu, double rho, double thickness,
			State ss, FluencyCriteria softeningRule, double alpha0, FunctionRnToR weightFunction,
			double interactionRadius, IntegrationPointIndexer indexer) {
		super(E, nu, rho, thickness, ss, softeningRule, alpha0);
		this.weightFunction = weightFunction;
		this.interactionRadius = interactionRadius;
		this.indexer = indexer;
	}

	public IsotropicDamageMMimplexNonlocal(IsotropicDamageMM mm, FunctionRnToR weightFunction,
			double interactionRadius, IntegrationPointIndexer indexer) {
		super(mm);
		this.weightFunction = weightFunction;
		this.interactionRadius = interactionRadius;
		this.indexer = indexer;
	}

	@Override
	public void Initialize(Element e) {
		super.Initialize(e);
		/*
		 * index the element
		 */
		indexer.index(e);
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		StructuralElement se = (StructuralElement) e;
		int dim = e.getMeshPart().getDimension();
		int ng = e.getNumIntegrationPoints();

		for (int k = 0; k < ng; k++) {
			double[] xi = e.getIntegrationPoint(k);
			IVector sigma = se.getStress(k);
			IVector corrected_sigma = getVariables(e).get(k).getVectorVariable(VariableName.SIGMA);
			IMatrix C = getConstitutiveMatrix(e, k);
			MutableDouble alpha = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA);
			double alpha_old = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA, -1);
			double alpha_old_old = getVariables(e).get(k).getScalarVariable(VariableName.ALPHA, -2);
			MutableDouble q = getVariables(e).get(k).getScalarVariable(
					VariableName.STRESS_LIKE_VARIABLE);
			double q_old = getVariables(e).get(k).getScalarVariable(
					VariableName.STRESS_LIKE_VARIABLE, -1);
			MutableDouble damage = getVariables(e).get(k).getScalarVariable(VariableName.DAMAGE);

			/**********************************************
			 * perform implicit update to compute corrected stress & damage
			 **********************************************/
			/*
			 * compute strain
			 */
			IVector epsilon = se.computeStrain(xi);

			/*
			 * compute effective stress
			 */
			BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, De, epsilon, 0.0, sigma);

			/*
			 * compute equivalent strain
			 */
			double etilde = 0;
			switch (equivStrainType) {
			case ENERGY_NORM:

				/*
				 * get list of neighbouring integration points
				 */
				double[] gpLoc = e.getMeshPart().XAt(xi);
				SpatialCriteria criteria = null;
				if (dim == 2) {
					criteria = new SphereSpatialCriteria(gpLoc[0], gpLoc[1], interactionRadius);
				} else if (dim == 3) {
					criteria = new SphereSpatialCriteria(gpLoc[0], gpLoc[1], gpLoc[2],
							interactionRadius);
				}

				if (integrationPointsMapping.get(e) == null) {
					integrationPointsMapping.put(e,
							new TreeMap<Integer, Map<Element, ArrayList<Integer>>>());
				}
				if (integrationPointsMapping.get(e).get(k) == null) {
					integrationPointsMapping.get(e).put(k, indexer.queryRange(criteria));
				}
				Map<Element, ArrayList<Integer>> ips = integrationPointsMapping.get(e).get(k);

				// if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				// Debugger.watch("checking neighbour element");
				// for (Element ie : ips.keySet()) {
				// Debugger.watch(ie.getId());
				// }
				// }

				/*
				 * iterate through each element and integration points lie
				 * within interaction domain
				 */
				double totalWeight = 0;
				for (Element ie : ips.keySet()) {
					ArrayList<Integer> points = ips.get(ie);
					for (int kk : points) {
						double[] tmpXi = ie.getIntegrationPoint(kk);
						double tmpW = ie.getIntegrationWeight(kk);
						double[][] jac = ie.getMeshPart().jacXAt(tmpXi);
						double detJ = LinAlgUtilities.determinantNxN(jac);

						/*
						 * compute equivalent strain at this point
						 */
						IVector tmpEpsilon = ((StructuralElement) ie).computeStrain(tmpXi);
						IVector tmpSigma = (IVector) tmpEpsilon.clone();
						BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, De, tmpEpsilon, 0.0, tmpSigma);
						double tmpEtilde = Math.sqrt(BLAM.dot(tmpEpsilon, tmpSigma));

						/*
						 * compute weight of nonlocal interaction
						 */
						double[] tmpGpLoc = ie.getMeshPart().XAt(tmpXi);
						double distance = LinAlg.diffNorm2(dim, gpLoc, tmpGpLoc);
						double nonlocalWeight = weightFunction.valueAt(distance);

						/*
						 * contribute to nonlocal equivalent strain and total
						 * weight
						 */
						etilde += tmpW * detJ * nonlocalWeight * tmpEtilde;
						totalWeight += tmpW * detJ * nonlocalWeight;
					}
				}

				/*
				 * scaling to achieve final nonlocal equivalent strain
				 */
				etilde /= totalWeight;

				if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
					Debugger.watch("totalWeight = ", totalWeight);
				}

				break;
			default:
				throw new Error("invalid equivalent strain type");
			}

			if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				Debugger.watch("etilde = ", etilde);
			}

			/*
			 * compute trial state
			 */
			double alpha_trial = alpha_old;

			/*
			 * check loading/unloading criteria (use g to be consistent with the
			 * paper)
			 */
			double g = etilde - alpha_trial;
			double dlambda;
			if (g <= 0.0) {
				/*
				 * unloading/elastic state
				 */
				dlambda = 0.0;
			} else {
				/*
				 * loading/damage state
				 */
				dlambda = g;
			}

			/*
			 * update internal variable
			 */
			alpha.setValue(alpha_old + dlambda);

			/*
			 * update stress-like variable
			 */
			double H = softeningRule.computeFluencyDerivative(alpha_old);
			// double H =
			// softeningRule.computeFluencyDerivative(alpha.getValue());
			q.setValue(q_old + H * dlambda);

			/*
			 * update damage variable
			 */
			damage.setValue(1 - q.getValue() / alpha.getValue());
			// if (damage.getValue() > 0.9999) {
			// damage.setValue(0.9999);
			// }
			if (damage.getValue() < 0.0) {
				/*
				 * just for debugging, to make sure the input & material model
				 * work correcly
				 */
				throw new Error("negative damage. It's unphysical");
			}

			/*
			 * update stress
			 */
			BLAM.copy(sigma, corrected_sigma);
			BLAM.scale(1 - damage.getValue(), corrected_sigma);

			/**********************************************
			 * perform explicit update to compute algorithmic stress and tangent
			 **********************************************/

			double dalpha = alpha_old - alpha_old_old;
			double alpha_alg = alpha_old + dalpha;

			// H = softeningRule.computeFluencyDerivative(alpha_old);
			double q_alg = q_old + H * dalpha;
			/*
			 * compute damage based on extrapolated value
			 */
			double damage_alg = 1 - q_alg / alpha_alg;
			/*
			 * compute extrapolated stress
			 */
			BLAM.scale(1 - damage_alg, sigma);
			/*
			 * compute extrapolated algorithmic tangent
			 */
			BLAM.copy(De, C);
			BLAM.scale(1 - damage_alg, C);

			/**********************************************
			 * debugging part
			 **********************************************/
			if (e.getEchoLevel().contain(EchoLevelType.DEBUG)) {
				Debugger.watch("finalize step:");
				Debugger.watch("epsilon =" + epsilon);
				Debugger.watch("sigma =" + sigma);
				Debugger.watch("alpha = ",
						getVariables(e).get(k).getScalarVariable(VariableName.ALPHA));
				Debugger.watch("damage = ",
						getVariables(e).get(k).getScalarVariable(VariableName.DAMAGE));
			}
		}
	}
}
